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QUALITATIVE AND ANALYTICAL RESULTS OF THE 
BIFURCATION THRESHOLDS TO HALO ORBITS 


SARA BUCCIARELLI, MARTA CECCARONI, ALESSANDRA CELLETTI, 
AND GIUSEPPE PUCACCO 

Abstract. We study the dynamics in the neighborhood of the collinear Lagrangian 
points in the spatial, circular, restricted three-body problem. We consider the case in 
which one of the primaries is a radiating body and the other is oblate (although the 
latter is a minor effect). Beside having an intrinsic mathematical interest, this model is 
particularly suited for the description of a mission of a spacecraft (e.g., a solar sail) to 
an asteroid. 

The aim of our study is to investigate the occurrence of bifurcations to halo orbits, 
which take place as the energy level is varied. The estimate of the bifurcation thresh¬ 
olds is performed by analytical and numerical methods: we find a remarkable agreement 
between the two approaches. As a side result, we also evaluate the influence of the differ¬ 
ent parameters, most notably the solar radiation pressure coefficient, on the dynamical 
behavior of the model. 

To perform the analytical and numerical computations, we start by implementing a 
center manifold reduction. Next, we estimate the bifurcation values using qualitative 
techniques (e.g. Poincare surfaces, frequency analysis, FLIs). Concerning the analytical 
approach, following [4j we implement a resonant normal form, we transform to suitable 
action-angle variables and we introduce a detuning parameter measuring the displace¬ 
ment from the synchronous resonance. The bifurcation thresholds are then determined 
as series expansions in the detuning. Three concrete examples are considered and we 
find in all cases a very good agreement between the analytical and numerical results. 
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1. Introduction 

We consider the motion of a small body with negligible mass in the gravitational field 
of two primaries which move on circular trajectories around their common barycenter. 
We refer to this model as the spatial, circular, restricted three-body problem (hereafter 
SCR3BP). As it is well known (see, e.g., [SU]), the SCR3BP admits five equilibrium 
positions in the synodic reference frame, which rotate with the angular velocity of the 
primaries. Two of such positions make an equilateral triangle with the primaries, while 
the other three equilibria are collinear with the primaries. 

While the triangular positions are shown to be stable for a wide range of the mass ratio 
of the primaries, the collinear points are unstable. Nevertheless, the collinear equilibria 
turn out to be very useful in low-energy space missions and particular attention has been 
given to the so-called halo orbits, which are periodic trajectories around the collinear 
points, generated when increasing the energy level as bifurcations from the so-called 
planar Lyapunov family of periodic orbits. 

We will consider a model in which one of the primaries (e.g. the Sun) is radiating; we 
will see that in some cases the effect of the solar radiation pressure (hereafter SRP) on 
a solar sail, namely an object with a high value of the area-to-mass ratio, is to lower the 
bifurcation threshold, enabling bifurcations before unfeasible. For sake of generality, we 
also consider the case in which the other primary is oblate, although this effect is definitely 
negligible in many concrete applications ([19]). We will consider the following three case 
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studies. The first one describes the interaction between the Earth-Moon barycenter and 
the Sun; we will refer to this case as the Sun-barycenter system. The second sample is 
provided by the Earth-Moon system. The last case describes the interaction between the 
Sun and one of the largest asteroids, Vesta, for which the effect of SRP is very important. 

Our model depends on three main parameters, which are the mass ratio /i of the 
primaries, the performance j3 of the sail providing the ratio between the acceleration of 
the radiation pressure to the gravitational acceleration of the main primary (equivalently, 
one can use the parameter q = 1 — (3) and the oblateness denoted by A, defined as 
A = J 2 r l, see [15, where J 2 is the so-called dynamical oblateness coefficient and r e is 
the equatorial radius of the planet. We remark that the orientation of the sail is set to 
be perpendicular with respect to the direction joining it with the main primary, so that 
the system still preserves the Hamiltonian character and no dissipation is allowed. All 
results could be easily generalised to the case in which the orientation of the sail is kept 
constant (i.e., non necessarily perpendicular) with respect to the direction joining it with 
the main primary as the system would still keep its Hamiltonian nature. 

The aim of this paper is to make concrete analytical estimates of the bifurcations to 
halo orbits based on an extension of the theory developed in [I] and to compare the 
mathematical results with a qualitative investigation of the dynamics. As a side result, 
we shall evaluate the role of the parameters of the model, in particular the solar radiation 
pressure coefficient and the oblateness. 

The collinear points with SRP and oblateness are shown to be of saddle x center x 
center type for typical parameter values (compare with Sections [2] and [3]). According to 
widespread techniques, the dynamics can be conveniently described after having applied 
a reduction to the center manifold, which allows us to remove the hyperbolic components. 
The procedure consists in applying suitable changes of variables to reduce the Hamilton¬ 
ian function to a simpler form and then to compute a Lie series transformation to get rid 
of the hyperbolic direction (see, e.g., mi, inn). This procedure is performed in Sections [3] 
and [5 taking care of the modifications required by the consideration of the SRP and of 
the oblateness. Once the center manifold Hamiltonian has been obtained, we proceed 
to implement some numerical methods to investigate the dynamics as the energy level is 
varied, precisely we compute some Poincare surfaces of section, we perform a frequency 
analysis and we determine the Fast Lyapunov Indicators anunusi). The independent 
and complementary application of the three methods has several advantages, as it helps 
to unveil many details of the complicated structures arising around the bifurcation values. 
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As for the analytical estimates (see 0), after the center manifold reduction we need 
also to construct a fourth order normal form around the synchronous resonance. Af¬ 
ter introducing a coordinate change to action-angle variables for the quadratic part of 
the Hamiltonian, we recognize the existence of a first integral of motion. Finally, we 
introduce a quantity called detuning, which measures the (small) discrepancy of the 
frequencies around the synchronous resonance. The bifurcation threshold will then be 
computed at the first order in the powers series expansion in the detuning (see Section [5]). 
We stress that, although we consider just the first non-trivial order of the expansion in 
the detuning, the analytical estimates are already in remarkable agreement with the nu¬ 
merical expectation (see Table |4j) . Clearly, more refined analytical results can be obtained 
computing higher order normal forms, at the expense of a bigger computational effort 
and provided the results are performed within the optimal order of normalization. 


To conclude, let us mention that our comparison is made on the three case studies 
mentioned before (Sun-barycenter, Earth-Moon, Sun-Vesta); however, we provide full 
details only for the Sun-Vesta case, since the other two samples have been already inves¬ 
tigated in the literature with an extensive use of the Poincare maps (see, e.g., BE]). 
The reason to focus on the case of Vesta is the following: due to the fact that Vesta has 
a relatively small mass, the interplay of the parameters gives rise to an interesting and 
non-trivial dynamical behavior (see Section 6.4). Indeed, in the Sun-Vesta case we shall 
see that the SRP plays an important effect, enabling more bifurcations at relatively low 
energy levels. 

In fact, in the purely gravitational model only the first bifurcation (the standard halo) 
occurs in general ([8]), when the frequency of the planar Lyapunov orbit is the same 
as the frequency of its vertical perturbation. A second bifurcation, in which the planar 
Lyapunov trajectory regains stability and a second unstable family appears, is a rather 
extreme phenomenon at high energy values. In the presence of SRP this second bifurca¬ 
tion occurs instead at a much lower value of the energy. Moreover, also a third bifurcation 
may occur in which the vertical Lyapunov loses stability and the second family disap¬ 
pears. We will show how it is possible to predict the influence of SRP on the thresholds 
with a first-order perturbation approach. 


This paper is organized as follows. In Section [2] we introduce the equations of motion 
and we determine the equilibrium points, whose linear stability as the parameters are 
varied is discussed in Section [3j where we prepare the Hamiltonian for the center manifold 
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reduction of Section [4} The analytical method to determine the bifurcation thresholds 
is presented in Section [5j while in Section [6] we implement the qualitative techniques, 
namely Poincare maps, frequency analysis and Fast Lyapunov Indicators. The compar¬ 
ison between the analytical and numerical approaches, as well as some conclusions, are 
given in Section [7} 


2. The three-body problem with an oblate primary and SRP 
In this section we start by introducing the equations of motion describing the SCR3BP 


with a radiating larger primary and an oblate smaller primary (see Section 2.1). Then, 
we proceed to compute the location of the collinear equilibrium points as a function of 


the parameters /i, (3, A (see Section 2.2). 


2.1. The equations of motion. Let us denote by /i and 1—/x the masses of the primaries 
with q G (0,1/2] and let us assume the units of measure, so that the gravitational 
constant is unity and the period of the primaries is equal to 2n. We consider a synodic 
reference system (O, X,Y, Z), rotating with the angular velocity of the primaries, with 
the origin located at the barycenter of the primaries and the abscissa along the primaries’ 
axis. The position of the smaller primary is at (—1 + /i, 0, 0), while the larger primary is 
located at (/x, 0,0). Let (X,Y,Z) be the coordinates of the third body in this reference 
frame. Let (Px, Py, Pz) be the conjugated kinetic momenta defined as Px — X — nY , 
Py — Y + nX, Pz = Z, where n denotes the mean motion, n = Zk jT rev and T rev is the 
period of revolution. 

The equations of motion are given by 


where 


ft = n(x, y, z) 


X - 2 nY 
Y + 2 nX 
Z 


<9ft 

~dX 

dY 

<9ft 

dZ ’ 


n 


y (X 2 + F 2 ) + 


q ( !-f) , p 


r i 


H- 

T2 




( 2 . 1 ) 


with the mean motion given by n — \J 1 + |A (compare with Appendix A), while 
the distances r 1 , r 2 from the primaries are given by r 1 = ■sj (A" — ji ) 2 + Y 2 + Z 2 and 
= A — /i + l) 2 + Y 2 + Z 2 (see, e.g., [5j HSJ). We have used q = 1 — (3, where the 
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LqQ 


performance /3 of the sail is defined by 

^ = 4t tc GM q B ’ 

where L 0 = 3.839 x 10 26 Watt is the Sun luminosity, Q is one plus the reflectivity, c is the 
speed of light, G is the gravitational constant, M 0 the Solar mass and B the mass/area 
ratio of the spacecraft. 

Notice that equations (2.1) are associated to the following Hamiltonian function: 

9(1 ~V) 


H(X,Y,Z,P x ,P y ,P z ) = -(Py + /y + P x ) + nY P x 


nXP, 


Y 


r i 


ri 


A 

1 2rf 


3Z 2 


( 2 . 2 ) 


2.2. The equilibrium points. It is well known (see, e.g., TO that the circular, re¬ 
stricted three-body problem admits the Lagrangian equilibrium points, precisely the 
triangular solutions, denoted as L 4 and L 5 , and the collinear solutions, denoted as Li, 
L 2 , L 3 , the latter being located on the axis joining the primaries. Due to the oblateness 
and to the radiation pressure effects, out-of-plane equilibria can be found in the spatial 
case (see 0). but we will not consider such solutions in the present work. 

To locate the equilibrium positions we impose that the derivatives of D in these points 
are zero, say |y = §y = §§ = 0. In particular, we obtain: 

qY( l-/i) 


d ft 2 V 

—— = n 1 
dY 


(y 2 + z 2 + (x — nYp 

3 AY fi 


+ 


15 AYZ 2 p 

2(Y 2 + Z 2 + (1 + X - p ) 2 ) 7 2 
Y/j, 


2 (y 2 + z 2 + { 1 + x - /i) 2 ) § (y 2 + z 2 + ( l + x- 

so that we have |p = 0 whenever Y = 0. In a similar way we obtain: 


p)- 


dn 

dz 


qZ(l~fi) 


(y 2 + z 2 + (x - /i) 2 )i 

9 AZp 


+ 


15 AZ 3 p 


2(Y 2 + Z 2 + (1 + X - p) 2 )^ 
Z fi 


2(Y 2 + Z 2 + { 1 + X - /i) 2 )5 (y 2 + z 2 + (l + X- /i) 2 )i 


and we have — 0 whenever Z — 0. As for the derivative with respect to the first 
component, we have: 

on 

dx 


= n 2 X 


q(l - p)(X - p) 


+ 


l5AZ 2 p{l + X — ji) 


(Y 2 + Z 2 + (X - p) 2 ) f 2(y 2 + Z 2 + (1 + X 
3A(l + X — p)p (l + X — p)p 




2 (y 2 + z 2 + ( l + x- p) 2 ) l (y 2 + z 2 + (i + x- p) 2 )l 


(2.3) 
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Inserting Y = Z = 0 in (2.3) we get the equation 

q{l — n)(X - n) 3An{l + X — n) 


n X — 


/u( 1 + X — /i) 


\X-^ 


= 0 . 


(2.4) 


2|1 + A" — yu| 5 |l + X-/i| 3 
Let us denote by 7 j the distance between Lj and the closer primary. Since the collinear 
points L 1 , L 2 , L 3 lie in the intervals (—1 + (— 00 , —1 + /r), (/i,+oo), setting X = 

7 i + h — 1 for Li, X = —72 + 11 — I for L 2 , X = 73 + h for L 3 , the quantity 7 j is found 
as the unique positive solution of the following generalized Euler’s equations: 

± 2n 1 2 7j + (2 n 2 /j, — 6n 2 )7® ± (6n 2 — 4n 2 /i)7| + (2 n 2 fi — 2 q/i — 2 n 2 + 2 q + 2 /i)^ 

+ 4 /rq^ =f (2/r + 3 A/x) 7 2 + 6^/17^ + 3 Afx = 0 


for j = 1 , 2 , where the upper sign holds for Li, while the lower sign holds for L 2 ; as for 
L 3) we have that 73 is the solution of the following Euler’s equation 


2 n 2 'jl + (8 n 2 + 2 n 2 fi)^ + ( 12 -n 2 + 8 n 2 /r )73 + ( 8 n 2 + 2 q/i + 12 n 2 fi — 2 q — 
+ (2 n 2 — 8 q — 4/i + 8 n 2 /i + 8 q/i )q| + (2 n 2 /! — 3 Afj, + 12 q/i — 12 q — 2^)73 
+ ( 8 q/i - 8^)73 + 2 qqi - 2 q — 0 . 


Making use of equation (2.4), we show below the dependence of the location of the 
collinear points as the parameters (/i, A, (3) are varied, where we assume that the param¬ 
eters belong to the following intervals 

0 < /r < 0.5 , 0 < A < 10 ~ 4 , 0 < 13 < 0.5 . (2.5) 


Here and in the following sections we consider three paradigmatic cases, which cor¬ 
respond to a spacecraft or a solar sail orbiting in the Earth-Moon system, in the Sun- 
barycenter system and in the Sun-Vesta system; these three cases encompass missions to 
satellites, planets or asteroids and are characterized by different values of the parameters 
as well as by different distances of the collinear points, as reported in Table [2j 


3. Linear stability of the collinear points and reduction of the 

QUADRATIC PART 

To study the stability of the collinear equilibrium points and their dependency on 
the parameters, we compute the linearization of the equations of motion, which requires 

1 The upper bound on A, say A < 10” 4 , is definitely large for solar system bodies, but it might apply 
to extrasolar planetary systems. 

2 A realistic upper bound on the sail performance should be (3 < 0.1; however, we consider f3 up to 0.5 
as it is often done in the literature, see e.g., [9]- 
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to expand the Hamiltonian in series np to the second order. Indeed, we will expand it 
directly up to the fourth order as this will be used for the reduction to the center manifold 
in Section |4} Although the computation of the linear stability is rather elementary, it is 
mandatory before performing the center manifold reduction as in Section [4j 


3.1. Expansion of the Hamiltonian to fourth order. We shift and scale the equi¬ 
librium points by making the following transformation (EH): 

X = T1jX + fi + a, Y = Tljy, Z = j j z, (3.1) 

where 7 j denotes again the distance of Lj from the closer primary; the upper sign cor¬ 
responds to L\ 2 and the lower sign to L 3 , while a = — 1 + 71 for L 1, a = —1 — 72 for 
L 2 and a = 73 for L 3 . These conventions will hold hereafter. It must be remarked that, 
in all cases, the transformation (3.1) is symplectic of parameter 7?, so that the resulting 
Hamiltonian must be divided by such factor. Inserting (3.1) in (2.3), we have: 

"g(l -/i) 


dVt , d 

dx = n + dx 


r 1 


h Afj, 
r 2 2 rf 


3y (z 2 A/j, 


2 f\ 


therefore, we obtain the equation of motion: 


± 2 ny ± n 2 x = 


n 2 (/i + a) Id 
^ 7 2 dx 


q( 1 ~ y) t h | Ay 
ri r 2 2 


37 ?z 2 A/i 


2 r\ 


di < 3 

where j = 1, 2, 3. We remark that the inverse of the distances are transformed as 
11 11 
r 1 


7 j 


{x =f ^) 2 + y 2 + - 2 


r 2 


7 j 


{—x ± g ^ L ) 2 + y 2 + z 2 


In a similar way, for the other two components we obtain: 


1 d 


=F7jZ/ =F 2n^j± = =F n (7 jy) T —~r 


li z = 


1 _ d_ 

7 jdz _ 


ij&y 

q( 1 - h) 


g(i ~h) + ii + 

r'i r 2 2r^ 

3rfiz 2 AfJl' 


37 fz 2 Afi 


2 r\ 


r 1 


h 71/i 
r 2 2r| 


2r| 


Expanding in Taylor series 1/ri, l/r 2 and their powers, we compute the linearized equa¬ 
tions as 


x — 2 ny — (n 2 — 2a)x = 0 
y + 2 nx + (— n 2 + b)y = 0 
z + cz = 0 , 

where the quantities a, b, c take different values according to the considered equilibrium 
point. Since in the following sections we shall analyze only L l and L 2 , we provide the 
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explicit values for such positions. Precisely, one has: 


a = 


b = 


c = 


q(l- H) 


cr 


+ 


h 


(1 + a)' 


T 


g(l ~ h) 

a 3 

q ( 1 - /0 


± 


± 




± 


3Ay 
(1 + a) 5 
3Ay 


(1 + a) 3 2(1 + a) 5 




± 


9A/j 


cr 


(1 + a) 3 2(1 +a) 5 ’ 


where the upper sign holds for Li and the lower sign for L 2 . 

Setting for shortness A = we can write a and c in terms of b as 

a — —(6 + A) , c — b + 2A . 

ft results that A > 0 and b > 0 for all A > 0, fi > 0, (3 > 0. 


(3,2) 


Finally, the complete equations of motion in the new variables can be written as 

— — H 

,,2 n 


x — 2ny — (n 2 — 2a)x 
y + 2 nx + (—n 2 + b)y 


Z + CZ = — 


7 2 <9x 

7? <%/ 
7 2 dz 


3 " n> 3 

1 

3 n >3 

1 

n> 3 


(3,3) 


where are suitable polynomials of degree n. Defining the momenta as p x = x — ny , 


p y = y + nx, p z = z and using (3.3), the Hamiltonian function (2.2) becomes 


H(x,y,z,p x ,p y ,p z ) = \{jp 2 x +p 2 y +p 2 z )+nyp x -nxp y +ax 2 + -by 2 + ) ) cz t ~\^ j H n . (3.4) 


7 


0 n >3 


3.2. Reduction of the quadratic terms and linear stability. The quadratic part 


of the Hamiltonian (3.4) is given by 

1 


H 2 {x,y,z,p x ,p y ,p z ) = ^{p 2 x + p 2 y ) + nyp x - nxp y + ^p 2 z + ax 2 + h)y 2 + ^cz 2 


By (3.2) it results that c > 0 for each of the equilibrium points. This implies that the 


vertical direction is described by a harmonic oscillator with frequency u> 2 = y/c, namely: 

dH 


Pz = 


dz 


= — CZ, 


dH 

Z = — = p z 

dp z 


As for the planar directions, following m the quadratic part of the Hamiltonian, that 
we keep denoting as H 2 , is given by 


1 ' „2 1 J 2 \ 1 „„„ 1 1 1 j 2 


H 2 (x, y,p x ,p y ) = ~(p x + p y ) + nyp x - nxp y + ax + -by- 
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Denoting by J the symplectic matrix, the equations of motion are given by 


/ x \ 

y 

Px 

\Py ) 

Let us define the matrix M as 


= JHess(H -2 


f x \ 

y 

Px 

\Py ) 


( o 


V 


n 1 
—n 0 0 

- 2 a 0 0 

0 —6 —n 


0 \ f x \ 

i y 

n p x 

0 J \ Py J 


M=J Hess{H 2 ) ; 

the associated characteristic polynomial is given by 

p( A) = A 4 + (2 n 2 + 2a + b)X 2 + (n 4 — 2 an 2 — bn 2 + 2 ab) 


(3.5) 


(3,6) 


Setting rj = A , the roots of the polynomial (3.6) are given by 


Vi = 


V2 = 


—2 n 2 — 2 a — b — a/16 an 2 + 4a 2 + 8 bn 2 — Aab + b 2 

2 

—2 n 2 — 2 a — b + a/16 an 2 + 4 a 2 + 8 bn 2 — 4 ab + b 2 


In order to study the stability of the collinear equilibrium points we have to establish 
the domains in which Ai ; 3 = and X 2 ,4 = ±y/V 2 are real, complex or imaginary. In 

particular, a given equilibrium point Lj will be linearly stable, if ip and p 2 are purely 
imaginary, while it is unstable elsewhere. We are interested to the case in which the 
collinear points are of the type saddlex center x center, which occurs whenever rp < 0 
and r ]2 > 0. This is equivalent to require that the following inequalities are satisfied: 

16an 2 + 4a 2 + 8 bn 2 — Aab + b 2 > 0 
—2 n 2 — 2a — b + a/16 an 2 + 4a 2 + 8 bn 2 — Aab + b 2 > 0 
—2 n 2 — 2 a — b — a/16 an 2 + 4a 2 + 8 bn 2 — Aab + b 2 < 0 . 


Using (3.2), the inequalities (3.7) become: 


9 b 2 + 4A(—4n 2 + A) + b(—8n 2 + 12A) > 0 
b - 2n 2 + 2A + a / 96 2 + 4A(-4n 2 + A) + b(-8n 2 + 12A) > 0 
b - 2n 2 + 2A - a/96 2 + 4A(-4n 2 + A) + 6(-8n 2 + 12A) < 0 . 


The first condition in (3.7) is satisfied whenever 


6 > -(2n 2 - 3A + 2y/n 4 + 6n 2 A) . 


(3.7) 


(3.8) 


Given that the function at the right hand side of (3.8) reaches its maximum n 2 for A = 


the inequality (3.8) is verihed if 6 > n 2 = \J 1 + | A. The second and third conditions in 
(3.7) can be reduced to verify just that 


86 2 + 6(8A — 4n 2 ) — 8n 2 A — 4n 4 > 0 , 
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which is satisfied again if b > n 2 , which holds for all values of the oblateness, the solar 
radiation pressure and the mass parameter in the intervals defined in (2.5). This con¬ 


cludes the discussion of the stability character of the collinear points, including the solar 
radiation pressure and the oblateness of one of the primaries. 


4. Center manifold reduction 

Due to the saddle x center x center character of the collinear equilibrium points (see 
Section [3]), we proceed to perform the reduction to the center manifold. The adopted 
procedure is a straightforward extension of that used in HU, provided that we include 
the necessary modifications to consider the effects of the oblateness and the solar radia¬ 
tion pressure. For completeness, we report here the main steps to treat the Hamiltonian 


(2.2) (see Appendix C for more details). We stress that after removing the hyperbolic 


direction, we will be able to perform a qualitative analysis based on Poincare sections, 
frequency analysis and Fast Lyapunov Indicators as it is done in Section [6] 

Taking into account that rji < 0 and r/ 2 > 0, let aq = \Z—pi an d Ai = we look for 
a change of variables, so that we reach a simpler form of the Hamiltonian. As described 
in detail in Appendix B, this is obtained by computing the eigenvalues and eigenvectors 


of M in (3.5), which provide a transformation allowing to get a Hamiltonian with the 


following quadratic part (with a slight abuse we keep the same notation for all variables): 

~ 2 1 ~ 2 ^ (4.1) 


H 2 (x,y,z,p x ,py,p z ) = \ixp x + —(py + y ) + + - ) , 


where 


Ai = 


UJl = 


-2 n 2 — 2 a — b + \/16 an 2 + 4a 2 + 8 bn 2 — Aab + b 2 


-2 n 2 — 2 a — b — a/16 an 2 + 4 a 2 + 8 bn 2 — Aab + b 2 


CU 2 = Vc • 

In analogy to HU, we introduce the following complex transformation 

x = qi, y = 

Px = Pi, Py = 


q-2 + W2 

V 2 ’ 

iq2 + P2 


z = 


V 2 


Pz = 


<?3 + m 

y /2 

'iq-s + P3 


V 2 
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which provides a complex expression for the Hamiltonian; we report here the quadratic 
part of the Hamiltonian, which takes the form: 


Hi(qi,q2,q3,PhP2,P3) = ^iqipi + i^iq2P2 + >: ~qq3P3 ■ 


(4.2) 


Beside the quadratic part (4.2), we need to compute the nonlinear terms. Straightforward 
but tedious computations, performed by means of the Mathematica® algebraic manipu¬ 
lator, allow us to find the expressions of H 3 and H 4 in complex variables. Afterwards we 
shall implement a Lie series transformation to obtain the reduction to the center manifold 
as described in the following section. 


4.1. Reduction to the center manifold. The reduction to the center manifold is 

obtained by making suitable changes of variables by using Lie series (see Appendix C). 
Indeed, we conjugate the SCR3BP to a Hamiltonian of the form 


H(q,p) = \1q1P1 + ioj 2 q 2 p 2 + ^ 3 ^ H n(<hP) > 

n> 3 

for suitable coordinates (gi, q 2 , q 3 ) and momenta (pi,P2,P3), where the quadratic part has 
been obtained in (4.2) and H n denotes a homogeneous polynomial of degree n. In the 
linear approximation the center manifold is obtained by imposing q\ = p\ = 0, since the 
hyperbolicity pertains to such variables. Then, we will require that g’i(0) = pi(0) = 0 
for gi(0) = pi(0) = 0, so that we obtain qi(t) = p\{t) = 0 for any time, due to the 
autonomous character of the problem. Taking into account Hamilton’s equations 


q% Hpi , Pi Hq, , 

this requirement is satisfied whenever in the series expansion of the Hamiltonian all 
monomials of the type hijq t pi with A 7^ ji are such that hij = 0, being i = (AA2A3) and 
j — jz)- In this way we obtain a Hamiltonian of the form H(q,p) = Hiy(q,p) + 

R N (q,p), where H N (q,p) is a polynomial of degree N in (q,p) without terms depending 
on the product qiPi, while R N (q,p) is a reminder of order N + 1. We refer to Appendix C 
for the description of a procedure based on Lie series to determine explicitly the required 
canonical transformation. Let us denote by (y, z,p y ,p z ) the normalized variables; the 
final expression of the Hamiltonian reduced to the center manifold has the following 
form: 

H(y, z,p yi p z ) = y (p 2 y + y 2 ) + y(p? + z 2 ) + H 3 {y,z,p y ,p z ) + H A {y,z,p y ,p z ) , (4.3) 
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where i/3, H 4 denote homogeneous polynomials of degree, respectively, 3 and 4. The 
frequencies cui, u; 2 , as well as those of the higher order terms, depend on the choice of 
the parameters and will be specified in each concrete case. 


5. Analytical estimates of the bifurcation thresholds 


Analytical results providing an estimate for the value of the thresholds at which the 
bifurcation of halo orbits takes place have been presented in [4]. The result is briefly sum¬ 
marized as follows. After the reduction to the center manifold, a normal form is computed 
around the synchronous resonance. The resulting normal form admits a first integral, 


related to the action variables of the harmonic oscillator (i.e., the quadratic part in (4.3)). 


A detuning measuring the displacement around the synchronous resonance is introduced; 
assuming that the detuning is small, one can compute the bifurcation threshold at differ¬ 
ent orders in the powers series expansion in the detuning. In [4J the computation at first 
and second order has been performed. Here, we extend the method of [4] by computing 
the thresholds for the model including oblateness and solar radiation pressure. As we will 
see in Section [6j the case in which the parameter /3 is different from zero allows one to 
find several bifurcations at relatively low energy levels. We anticipate that the analytical 
estimates computed in this section will agree with the numerical values of Section [6] (see 
Table g. 


Let us write (4.3) in complex form, implementing the following change of coordinates: 

y/2 V2 

Q 2 = ~^r(P2+iq2) , Qs = -^-(P3 + * 93 ) , 


P2 = -K-I2) , 


Pi = - iqs) 


Next, we introduce action-angle variables (i 2 ,i3, #2, #3) for the quadratic part of the 
Hamiltonian by means of the change of coordinates: 


Q 2 = -iyfl,~oe id2 


Q 3 = -iy/h e i9s 


o —i02 


p 3 = y/h e ~ id3 


P2 — \fl 2 e 

so that we obtain a Hamiltonian of the form 

H( 02 , 0'3i h, I 3 ) ~ U 2 I 2 + U 3 I 3 + H 3 ( 02 , O 3 , i 2 , I,) + H i { 92 ^ 63 , i 2 , i 3 ) • 


To investigate the appearance of the resonant periodic orbits we compute a normal form 
in the neighborhood of the synchronous resonance cu 2 = 003 . The explicit computation 
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TABLE 1. Coefficients of the normal form (5.1) around L x for the Earth- 
Moon, Sun-barycenter, Sun-Vesta systems with and without solar radia¬ 
tion pressure. 



020 

002 

Oil 

bn 

Earth-Moon 

0.162109 

0.144891 

0.0726274 

0.116537 

Sun-barycenter (3 = 0 

0.0989667 

0.08098 

-0.0256235 

0.102138 

Sun-barycenter f3 = ICE 2 

0.136123 

0.108011 

0.0189407 

0.11106 

Sun-Vest a (3 = 0 

0.0957347 

0.0777063 

-0.0304854 

0.101319 

Sun-Vesta (3 = HE 2 

0.0157472 

0.00203253 

4.11966 10“ 7 

0.00533371 


shows that the Erst non-trivial order is given by H±, since the third degree term // 3 does 
not contain resonant terms. We are thus led to a resonant normal form given by 


Hnf(@2, @ 3 i h, I 3 ) — ^ 2 I 2 + tU 3 / 3 + 0 2 q/ 2 + a 02^3 + -^-^(Oll + 2 &n cos(20 2 — 2# 3 ) , (5.1) 


where the coefficients a 2 0, o 02 , an, bn are evaluated explicitly in Table [l] 

The dynamics is determined by the normal modes Ik = const., k — 2,3 and by the 
periodic orbits in general position, related to the resonance. The normal modes always 
exist and at low energies are both stable: I 2 = const., J 3 = 0 gives rise to the planar 
Lyapunov orbit; J 2 = 0, / 3 = const, gives rise to the vertical Lyapunov orbit. When 
stable, these orbits are surrounded by families of Lissajous tori. The resonant families 
may appear as bifurcation from the normal mode at some given energy threshold and 
are determined by the condition that the frequency of the normal mode is equal to that 
of its normal perturbation. The normal modes can be again stable through a second 
bifurcation; a concrete example of a second bifurcation for the Earth-Moon case is given 
in [8] (we refer to pfl UHlUTl 18] for further details). To investigate the possible sequences 
of bifurcations we proceed as follows. 


From Hamilton’s equations associated to (5.1), it is readily seen that J 2 + J 3 = 0. This 
remark leads to introduce the following change of variables ([4]): 


£ — I'2 + 1.3 , K-h, 

» = 0 3 , ^ = e 2 -e 3 . (5.2) 

Moreover, following (3J let us introduce the detuning 6 as 

5 = uj 2 - cu 3 , 
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which measures the displacement from the synchronous resonance. Using (5.2) the Hamil¬ 
tonian (5.1) becomes 

H ne w(£, H, v, VO = £ + 81Z + alZ 2 + b£ 2 + cElZ + d(TZ 2 — EH) cos 2ijj , 

where 8 = 8 /M 3 , a = (gqo + a02 — On)M, 5 = CI 02 /U 3 , c = (an — 2ao2)/<V3, d = — 2 bn/a> 3 - 
It has been shown in [4] that the equilibria associated to Hamilton’s equations of H new 
can be classified as inclined (or ‘anti-halo’) when if} = 0 or ip = tt, and loop (or ‘halo’), 
when -0 = ±7 t/ 2. In view of the reflection symmetries, each case actually corresponds to 
a double family. They exist at the following level values of the integral of motion E: 


E, = Sujl 

—an + 2(020 — bn) 


for the inclined families and 

E(,y = 


&4 


—an + 2(a2o + &11) 


£. = 


Eiz — 


8cu?, 


—2oq2 + an + 2&n 


—2a Q 2 — On + 2bn 


(5.3) 


(5.4) 


for the loop families. These values correspond to a first order computation in the detun¬ 
ing; refined (but more complicated) values at second order have been found analytically 

in [4]. 

The physical interpretation of the thresholds in (5.3), (5.4) is the following. The 
quantity Ee y determines the bifurcation of the halo families from the planar Lyapunov 
orbit, when this becomes unstable. At Ei y the planar Lyapunov orbit turns back to 
being stable and one observes the bifurcation of the (unstable) anti-halo orbits. Finally, 
the two unstable families which have been formed collapse on the vertical at E lz . The 
disappearance of the halo at Ei z never occurs in all cases we have investigated. 

This sequence of bifurcations will be clearly shown in the case of the asteroid Vesta 
under the effect of solar radiation pressure (see Section [6]). 


6. Qualitative analysis of the bifurcation values 

On the basis of the center manifold reduction obtained in Section [4| we implement some 
numerical techniques which allow us to characterize the dynamics and, in particular, to 
distinguish between the different types of orbits, precisely planar Lyapunov and halo 
orbits, the latter ones obtained at specific levels of the energy at which the bifurcation 
takes place. As concrete models, we consider three paradigmatic cases, characterized by 
different values of the mass ratio: a relatively high value as in the Earth-Moon system, 
an intermediate mass ratio as for the Sun-barycenter system (between the barycenter of 
the Earth-Moon system and the Sun), and a low value as in the Sun-Vesta system. 
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6.1. Poincare section. To get a qualitative description of the dynamics in the center 
manifold we start by computing a Poincare section as follows. We set z — 0 and we 
fix an energy level H = ho, from which we compute the initial value of p z choosing the 
solution with p z > 0. The Poincare section is then shown in the plane ( y,p y ); we will see 
that, as the energy increases and exceeds a specific energy value, halo orbits arise from 
bifurcations of planar Lyapunov periodic orbits. 


6.2. Frequency analysis. This technique consists in studying the behavior of the fre¬ 
quency map mm), which is obtained computing the variation of the absolute value 
of the ratio of the frequencies, say u r = \u y /u z \, as a function of the initial values of the 
action variables, whereas the initial conditions of the angles can be set to zero (see, e.g., 
[3], see also [2]). 

The frequency analysis has the advantage to be computationally fast and it allows us 
to obtain a complementary investigation of the occurrence of halo orbits. Precisely, we 
proceed as follows. Concerning the initial conditions, we fix as starting values z — 0 
and p y = Py, we scan over initial values for y in a given interval and for an assigned 
energy level H = ho, we compute the corresponding value of p z . We find convenient to 
avoid using Cartesian variables, and we rather transform to action-angle variables for the 
quadratic part of (4.3). Thus, we introduce harmonic oscillator actions ( J y ,J z ) defined 
through the expressions 


p y = y/2Jy cos 9 y , y = \plJy sin 9 y , 

Pz = V2J~z cos 9 Z , z = \/2T z sin 9 Z , 


where {9 y , 9 Z ) denote the conjugated angle variables. Next we perform a first order pertur¬ 
bation theory by averaging over the angle variables to obtain a normalized Hamiltonian, 
whose derivative provides an expression for the frequencies associated to the given ini¬ 
tial conditions. Finally, we back-transform to the variables (J y , J z ) to get a frequency 
vector ( in y ,uj z ) associated to the previous initial data. The frequency map is obtained by 
computing the variation of u r = \uj y /uj z \ as the initial condition is varied. 


6.3. Fast Lyapunov Indicator. In order to investigate the stability of the dynamics in 
the center manifold, we compute a quantity called the Fast Lyapunov Indicator (hereafter 
FLI), which is determined as the value of the largest Lyapunov characteristic exponent 
at a fixed time (see m By comparing the values of the FLIs as the initial conditions 
or suitable parameters are varied, one obtains an indication of the dynamical character 
of the orbits (precisely Lyapunov or halo) as well as of their stability. The explicit 
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TABLE 2. Main parameters and location (in normalized units) of the 
collinear points for the Earth-Moon, Sun-barycenter, Sun-Vesta systems. 



Earth-Moon 

Sun-barycenter 

Sun-Vest a 

p 

1.2154 • 10" 2 

3.040423 ■ 10~ 6 

1.3574 • 10" 10 

J 2 

2.034- 10“ 4 

1081 • 10" 6 

0.0812232 

A 

4.15559 • 10~ 9 

1.96782 • 10- 12 

4.54776 • 10" 14 

P 

0 

10- 2 

O 

1 

to 

L\ 

-0.836898 

-0.988731 

-0.996651 

L 2 

-1.1557 

-1.00908 

-1.000115 

l 3 

1.00506 

0.996657 

0.996655 


computation of the FLI proceeds as follows. Let £ = (y, z,p y ,p z ), let the vector field 


associated to (4.3) be denoted as 


£ = /(£), £ e R 4 , 

and let the corresponding variational equations be 

!Z= r l ’ !Z eR • 

Having fixed an initial condition £(0) G R 4 , ?/(0) G R 4 , the FLI at a given time T > 0 is 
obtained by the expression 

FLI(£(0), 77 ( 0 ), T) = sup log ||/ 7 (t)|| , 

0 <t<T 

where || • || denotes the Euclidean norm. 


6.4. Applications. We proceed to implement the techniques described in Section 6.1 


6.2, 6.3 to the concrete samples provided by the Earth-Moon, Sun-barycenter, Sun- 


Vesta systems. The computations have been performed using Mathematica® as well 
as developing dedicated programs in a general-purpose programming language. The 
parameters associated to these three samples are listed in Table [2] 

The values of the quantities introduced in Section [4| needed for the center manifold 
reduction, are listed in Table [3j 


We analyze in detail the Sun-Vesta case, which presents several interesting features 
as the different parameters are varied. We start by computing the Poincare surfaces 
of section of the center manifold associated to L^. We report in Figure [T| the Poincare 





















18 


S. BUCCIARELLI, M. CECCARONI, A. CELLETTI, AND G. PUCACCO 


Table 3. Quantities for the center manifold reduction associated to the 
Earth-Moon, Sun-barycenter, Sun-Vesta systems. 



Earth-Moon 

Sun-barycenter 

Sun-Vesta 

7i 

0.150948 

0.01127 

0.00334854 

«( 7i) 

-0.836898 

-0.98873 

-0.996651 

a 

-5.14772 

-3.15056 

-1.00363 

b 

5.14772 

3.15056 

1.00363 

c 

5.14772 

3.15056 

1.00363 

Ai 

2.9321 

2.13994 

0.10407 


2.33441 

1.85169 

1.0036 

LU2 

2.26886 

1.77498 

1.00181 

Si 

14.9084 

9.6584 

0.79682 

S 2 

23.4324 

12.6138 

2.02523 


sections in the plane (■ y,p y ) with only the gravitational effect, namely with j3 — 0 (no 
solar radiation pressure) and A — 0 (no oblateness). The maps show that the appearance 
of halo orbits occurs for an energy value approximately equal to h = 0.35 (more accurate 
computations will provide the bifurcation vale h = 0.3341, compare with Table [4]). For 
higher levels of the energy, the amplitude of the halo orbits increases as shown by the 
map at h — 0.5. 

A different situation occurs when the radiation pressure and the oblateness are switched 
on, as shown in Figure [2] which reports the Poincare sections for the case with (3 = 10~ 2 
and A = 4.54776 • 10 -14 . Indeed, we have noticed that the oblateness has a small effect, 
while the parameter /3 plays a major role in shaping the dynamics. In fact, a comparison 
between Figure [T| and Figure [2] shows that the sequence of bifurcations is completely 
different. 

In Figure [2] we observe a regular behavior for h = 0.04, while already at h = 0.05 the 
dynamics experiences a first bifurcation with the appearance of halo orbits and simul¬ 
taneous loss of stability of the planar Lyapunov orbit. A second bifurcation of unstable 
inclined orbits takes place at about h = 0.1; at this stage, the planar Lyapunov orbit 
regains stability (compare also with | 8 j), as shown in Figure [2] where the planar Lyapunov 
orbit is given by the outermost curve. For increasing values of the energy, the unstable 
families (which are located on the vertical axis of the plot for h = 0.4) collapse on the 
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Figure 1. Poincare sections of the center manifold associated to Li of 
the Sun-Vesta system on the plane ( y,p y ) without solar radiation pressure 
(i.e., (3 = 0) and no oblateness (i.e., A = 0); different values of the energy 
are taken into account: h = 0.2 (left panel), 0.35 (middle), 0.5 (right). 



Figure 2. Poincare sections of the center manifold associated to L\ of the 
Sun-Vesta system on the plane ( y,p y ) with /3 = 10~ 2 and A = 4.54776 ■ 

10~ 14 ; different values of the energy are taken into account: h = 0.04 
(upper left panel), 0.05 (upper right), 0.1 (bottom left), 0.4 (bottom right). 

vertical Lyapunov orbit at the center of the axes (see the bottom right panel of Figure |2|. 
This corresponds to the third bifurcation, which is shown at about h = 0.4. 

The above results are confirmed by the study of the model through frequency analysis, 
as shown in Figures [3] and [4| 

In particular, in Figure [3] we provide the results of the frequency analysis for the Sun- 
Vesta case without solar radiation pressure and no oblateness in the (J°,oy) plane with 
Jy the initial condition and oj r = \u} y /uj z \. The first plot corresponds to h = 0.2 and it 
shows a regular behavior as it was found in the first plot of Figure [lj The occurrence 
of small halo orbits in the frequency analysis investigation corresponds to the two tiny 
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/ 0.0104 \ 

i/o.4 
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' v 0.0100 


\ Jy 0 
0.2 0.4 ■OA* 


Figure 3. Frequency analysis in the plane (J°, ay) for the Sun-Vesta case 
with (3 = 0, A = 0; different values of the energy are taken into account: 
h = 0.2 (left panel), 0.35 (middle), 0.5 (right). 


bumps at the outermost sides of the plot for h = 0.35 in Figure [3j these bumps increase 
in size for h = 0.5, in full agreement with the corresponding plot of Figure [TJ 

The results of the investigation through frequency analysis in the case with solar ra¬ 
diation pressure and oblateness are provided in Figure |4[ which corresponds to the Sun- 
Vesta case with (3 = 1CT 2 and A = 4.54776 • 10~ 14 . Again we find full agreement with the 
Poincare maps provided in Figure [2} Precisely, the upper left panel of Figure [4] shows a 
regular behavior, while tiny bumps, corresponding to the bifurcation of the halo orbits 
for h = 0.05, are present in the upper right panel of Figure |4j The three island regimes 
occurring for h = 0.1 correspond to the central bump and the two left and right wings of 
the bottom left panel. Finally, for h = 0.4 we observe a singular behavior on the vertical 
axis of the bottom right panel of Figure [4], which corresponds to the vertical Lyapunov 
orbit at the origin of the coordinates in Figure [2] 

The computation of the FLIs provides additional information: beside the overall struc¬ 
ture of the phase space, it yields the regular or chaotic character of the different trajec¬ 
tories. In particular, we can locate the separatrices, which were not easily determined 
within the Poincare maps (compare, e.g., Figure [2] bottom left and Figure [6] bottom left). 

In Figure [5] we provide the results obtained computing the FLIs for the Sun-Vesta 
case without solar radiation pressure and no oblateness; the results must be compared 
with Figure [T| in order to distinguish the different orbits on the Poincare maps and to 
determine their stability on the FLI plots. The color bar on the side of each plot gives 
the quantitative value of the FLI. 

In Figure |6]we provide the results obtained computing the FLIs for the Sun-Vesta case 
with (3 = 10~ 2 and A = 4.54776-10 -14 . Again, Figure[6]must be compared with Figure[2] 
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Figure 4. Frequency analysis in the plane (J°, uj y ) for the Sun-Vesta case 
with (3 = 10 -2 and A = 4.54776 • 10~ 14 ; different values of the energy are 
taken into account: h = 0.04 (upper left panel), 0.05 (upper right), 0.1 
(bottom left), 0.4 (bottom right). 



Figure 5. Fast Lyapunov Indicators for the Sun-Vesta case with f3 = 0, 
A = 0; different values of the energy are taken into account: h = 0.2 (left 
panel), 0.35 (middle), 0.5 (right). 


in order to identify the various trajectories on the Poincare sections and to characterize 
their stability on the FLI plots. It is remarkable how the FLI plots highlight also the 
separatrices as shown in particular in the bottom panels of Figure [6] 

When the halo orbits are very tiny, a zoom becomes necessary, as shown in Figure [7J 
where we provide a magnification of the plots obtained in the cases with /3 — 0, h — 0.35 
(left panel) and f3 = 10~ 2 , h = 0.05 (right panel). 

In the panels of Figure [6] we notice some lighter regions which are mainly along the 
horizontal direction (compare with the bottom right panel); these zones do not have a 
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Figure 6. Fast Lyapunov Indicators for the Sun-Vesta case with /3 = 1CF 2 
and A = 4.54776 • 10 -14 ; different values of the energy are taken into 
account: h = 0.04 (upper left panel), 0.05 (upper right), 0.1 (bottom left), 
0.4 (bottom right). 
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Figure 7. A zoom on the cases with /3 = 0, h = 0.35 (left panel) and 
/3 = 10” 2 , h = 0.05 (right). 

real dynamical meaning, but they are rather an artefact of the choice of the initial tangent 
vector, as the FLI strongly depends on this choice. In Figure [6] the tangent vector has 
been fixed as (v Px ,v Py ,v x ,v y ) = (1,0, 0,0) and we observe that lighter regions occur in 
the direction perpendicular to the chosen tangent vector. A reliable description of the 
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TABLE 4. Numerical (num) and analytical (anal) value of the energy at 
which the bifurcation to halo orbits has taken place; the results are given 
for the Earth-Moon, Sun-barycenter, Sun-Vesta systems. 



L | num 

L\ anal 

L ) num 

L 2 anal 

Earth-Moon f5 = 0, A = 0 

0.3026 

0.3069 

0.3776 

0.3636 

Earth-Moon j3 = 0, A = 4.15559 ■ 10 -9 

0.3026 

0.3069 

0.3776 

0.3636 

Sun-barycenter f3 = 0, A = 0 

0.3333 

0.3356 

0.3377 

0.3391 

Sun-barycenter /3 = 1CT 2 , A = 1.96782 • 1CT 12 

0.2793 

0.2864 

0.3759 

0.3755 

Sun-Vesta /3 = 0, A = 0 

0.3341 

0.3373 

0.3346 

0.3374 

Sun-Vesta f3 = 1CT 2 , A = 4.54776 • 1CT 14 

0.0422 

0.0424 

0.4451 

0.4434 


dynamical character of a specific region by means of the FLIs can only be obtained by 
comparing the plots produced using different tangent vectors in orthogonal directions or, 
alternatively, by increasing very much the accuracy of the computations. However, we 
remark that the analysis of Figure [6] suffices to distinguish the bifurcations as well as the 
main structures, like halo orbits or separatrices; further refinements go beyond the aims 
of the present work. 


7. Analytical versus numerical results 


In this section we compare the results which are obtained implementing the analytical 


formulae (5.3)-(5.4) for the computation of the bifurcation thresholds with the numerical 


values obtained through the Poincare maps or the FLIs. 

The analytical and numerical bifurcation values for all three case studies (Earth-Moon, 
Sun-bary center, Sun-Vesta) are listed in Table [4j The analytical results require the 
computation of the normal form as well as the computation of the thresholds at first 


order in the detuning as given by (5.3)-(5.4). The qualitative results are obtained looking 


at the bifurcations observed on the Poincare sections as well as through the FLI maps 
(the results are also validated through the application of the frequency analysis). 


From Table [4] we can draw the following conclusions. 

(i) The agreement between the analytical and numerical results is very satisfactory. 
The relative error between the theoretical and experimental values ranges between 
10~ 2 and 10 -3 . 



















24 


S. BUCCIARELLI, M. CECCARONI, A. CELLETTI, AND G. PUCACCO 


(ii) A first order estimate is already enough to get a good approximation of the 
bifurcation thresholds; this estimate requires a very little computational effort 
with respect to the qualitative analysis based on the Poincare maps or the FLIs. 
Obviously, better results can be obtained computing higher order normal forms, 
but at the expense of dealing with more complex formulae. 

(in) Switching on the solar radiation pressure provokes drastic changes for small mass 
parameters. In particular, in the Sun-Vesta case the first, second and third 
bifurcations take place at much lower values of the energy level, such that the 
other bifurcations become feasible. From the physical point of view, the reason 
for such a peculiar behavior is due to the balance between a smaller mass like 
that of an asteroid and the effect of the solar radiation pressure. 

(iv) The role of the oblateness is essentially negligible in all considered cases. This fact 
could have been inferred easily, but we believe worthwhile to derive a complete 
model, valid not only for the cases studied in the present paper, but also for 
general situations in which the small body could have a very irregular shape. 
Simple experiments show that in a Sun-asteroid sample, the oblateness becomes 
important only when the factor A is as large as 10 -6 — 10 -7 . This parameter 
value does not apply to Vesta, but it might be of interest for other astronomical 
situations. 
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Appendix A. Derivation of the mean motion around an oblate primary 

For completeness we derive here the formula for the mean motion of a massless body 
P under the gravitational effect of an oblate primary. Let us denote by Mp and A the 
mass and the oblateness coefficient of the primary, r represents the distance of P from 
the primary and P is the angular momentum constant. Then, the effective potential 
(see, e.g., HI) can be written as 

n 2 M P M P A 

2 ^ 


V'„(r) 


2 r 3 
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whose derivative is 


V eff 


, , U 2 Mp 3 MpA 

j) =-T + —5“ + 


2 r 4 


The solutions of K//( r ) = 0 are given by 


r = 


U 2 ± - 6M 2 A 


2 Mr 


(A.l) 


Taking into account that A is small, we can approximate the non-trivial solution of (A.l) 
by 

PL 2 3 Ml A, 


r = 




(A.2) 


m p 2 n 2 ' 

Assume that the orbit of P is circular, say r = a, we have that PL 2 = n 2 a 4 which, 


together with (A.2) and the normalization of the units of measure such that a — 1 , 


Mp = 1, provides: 


n 2 = 1 H —A . 


Appendix B. Reduction of the quadratic part 
In this section we provide the details of the reduction of the quadratic part as per¬ 


formed in Section 3.2 in order to obtain the Hamiltonian (4.1) (equivalently (4.2)). The 
procedure is very similar to that explained in HU, to which we refer for a complete discus¬ 
sion; however, for self-consistency, we provide here some details containing the necessary 
amendments to encompass the oblate case with solar radiation pressure. 


We start by computing the eigenvalues of M in (3.5); we denote by I n the n-dimensional 


identity matrix. We notice that by defining M\ = M — A/4, we can write M\ = 
^ A ) f ^ n \ J and B = ( ^ a j. Then, the kernel of 


0 -b 

Wl 

w 2 

computations show that the eigenvector of M is given by 


A J \ —n -A J \ 

M\ is given by the solution of M\w = 0 with w = 


and ivi,w 2 G R 2 . Simple 


(2nA, A 2 + 2a — n 2 , n\ 2 — 2an + n 3 , A 3 + (2a + n 2 )A) T , 
where A is an eigenvalue of M (the superscript T denotes the transposed). 


Let us consider the eigenvectors associated to oj\ = y/—rj 1; from (3.6) we have p( A) = 0, 
so that uj\ satisfies the equation 


cu 4 — (2n 2 + 2a + b)u\ + (n 4 — 2 an 2 — bn 2 + 2 ab) = 0 
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Using iui = Ai, the eigenvector u Ul + iv Wl associated to is given by 
u Wl + iv Ul = (2niu>i, —uf + 2 a — n 2 , 

— mol ~ % an + n 3 , —iu\ + (2a + n 2 )iui) T , 


while for ±Ai = we obtain: 

u + y = (2nAi, A 2 + 2a — n 2 , 

n\\ — 2 an + n 3 , Xl + (2a + n 2 )Ai) T 
V-\ l = (—2nAi, A 2 + 2a — n 2 , 

nA 2 — 2 an + n 3 , —A^ — (2a + n 2 )Ai) T . 


Let C = {u +Xl , v» U-Ai, Uwi); we have 

C T JC = ^ _° D ^ , where D 
with d Xl and d m given by 


d Xl 0 A 

0 d u}1 J 


d Xl = —2Ai((—4 n 2 + 2 a — 6)A 2 — 4n 4 + bn 2 + 6 an 2 — 2 ab + 4a 2 ) , 
duj i = —cai((—4n 2 + 2a — 6)ca 2 + 4n 4 — bn 2 — 6cm 2 + 2 ab — 4a 2 ) . 

In order to obtain a symplectic change of variables, we re-scale by Si = \Jd Xl , s 2 = \/d ui 
and we require that d Xl > 0, d Wl > 0. Finally, we re-scale (z,p z ) by (-T=, ^/Eq)- The 
hnal symplectic change of variables is given by the matrix whose columns are u +Xl /s i, 
U u>l / s 2i V_ Xl /si, Vlji /52- 


Appendix C. Center manifold reduction 


Let H be a Hamiltonian function with 3 degrees of freedom, admitting an equilibrium 
point of type saddlexcenterxcenter. Let ±A, ±iuq, ±iw 2 be the eigenvalues of the 
linearized system. Expanding the Hamiltonian around the equilibrium point in complex 
variables, we obtain a simpler Hamiltonian function of the form 


where 


H(q,p) = H 2 (q,p) + ^H n (q,p) , 

n> 3 


H 2 {q,p) = Xqipi + iu)iq 2 p 2 + iu 2 q 3 p 3 , 

while H n (q,p ) are homogeneous polynomials of degree n in the variables (q,p). To decou¬ 
ple the hyperbolic direction from the elliptic one, we need to kill the monomials whose 
exponent in p\ is different from that in q\. This procedure, which makes use of Lie series, 
allows one to get a first integral with level surface given by the center manifold. We 
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sketch below the procedure to find the canonical transformation, referring to [11] for full 
details. 


Given a Hamiltonian H and a generating function G, we denote by H the functioij^] 
H = H + {H, G} + i {{H, G} , G} + | {{{ff, G} , G} , G} + ... (C.l) 

If G has degree 3, say G = G 3 . comparing same orders in ( |C.l ) provides 

H 2 = h 2 , 

H, = H 3 + {H 2 l G 3 } , 

Hi = H 4 + {H 2 , G 4 } + {H 3 ,G 3 } + — {{H 2 ,G 3 } ,G 3 } , 

Next we look for G 3 such that H 3 is in normal form. Expanding H 2 , H 3l G 3 as 

3 

H 2 (q,p) = ^VjQjPj , 

3 = 1 

H 3 (q,p) = h k q ,k p Q kq P kp , 

| kq | | kp | = 3 

G 3 {q,p) = 9k q ,k p q kq p kp 

| kq | -hi Aip |=3 

for some coefficients hk q ,k v , 9 k q ,k p , Vj and denoting by 7 ] = (\,iui,iu 2 ), we have 


-h 


'knMn h h 

q p q q p p 


G 3 (q,p)= J2 < k -k v> 

(k q ,k p )es 3 K kp ^ 11 > 

where S 3 is the set of indexes (k p , k q ), such that |/c p | + |A: g | =3 and with the first component 

of k p different from the first component of k q . We proceed iteratively to higher orders up 

to a given order, say N, so that we obtain: 

H(q,p) = H (N \q 1 p 1 ,q 2 ,q 3 ,p 2 ,p 3 ) + R (N+1 \qi ,q2,q3,pi,P2,p 3 ) , 


where is a polynomial of degree N and f?^ JV+1 ^ ) is a reminder of order IV + 1. Ne¬ 
glecting the reminder and setting q\p\ = 0, we eliminate the hyperbolic component of 
H^ N > and we obtain a 2 degrees of freedom Hamiltonian of the desired form. As remarked 
in mi, there are no small divisors in the above procedure, since | < k q — k p , v > \ > |Ai| 
for any (k p , k q ) e Sj, j > 3. 


3 Curly brackets denote, as usual, the Poisson brackets ( 0 )- 
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Table 5. Coefficients up to degree 4 of the Hamiltonian restricted to the 
center manifold for the Sun-Vesta system with f3 = 1CT 2 . The exponents 
(ki, k 2 , ks, k 4 ) refer to the variables (y, z,p y ,p z ). 


k i 

k 2 

h 

k 4 

h k 

2 

0 

0 

0 

0.501797549378742 

0 

2 

0 

0 

0.500906031584819 

0 

0 

2 

0 

0.501797549378742 

0 

0 

0 

2 

0.500906031584819 

2 

0 

1 

0 

0.0014920550494420622 

0 

0 

3 

0 

-0.000248666270046148 

0 

2 

1 

0 

0.0003790464810461912 

4 

0 

0 

0 

-0.02099512477285749 

2 

2 

0 

0 

-0.010667382077111268 

0 

4 

0 

0 

-0.001354993618855492 

2 

0 

2 

0 

0.04199066782897781 

0 

2 

2 

0 

0.010667253491139606 

0 

0 

4 

0 

-0.003498979471471415 

1 

1 

1 

1 

2.81508155360122 • 10~ 7 

2 

0 

0 

2 

1.8824017340799554 • 10" 7 

0 

2 

0 

2 

4.7821141283422436 • 10" 8 

0 

0 

2 

2 

-9.411646402154861 ■ 10~ s 


Going back to real variables {y, z,p y ,p z ), we obtain a Hamiltonian of the form 
H(y, z,p y ,p z ) = Y hkiMfoM y kl z k2 p k yP k z ■ 

k\ ,/C2,fc3,fc4GZ 

The hrst few non-zero terms h k] .k 2 M-M °f the Hamiltonian restricted to the center man¬ 
ifold are provided in Table [5j 


References 

[1] A. Celletti, Stability and Chaos in Celestial Mechanics , Springer-Verlag, Berlin; published in associ¬ 
ation with Praxis Publishing Ltd., Chichester, ISBN: 978-3-540-85145-5 (2010) 

[2] A. Celletti, C. Froeschle, On the determination of the stochasticity threshold of invariant curves, Int. 
J. Bif. Chaos 5 , n. 6, 1713-1719 (1995) 

[3] A. Celletti, C. Froeschle, E. Lega, Frequency analysis of the stability of asteroids in the framework 
of the restricted, three-body problem, Cel. Mech. Dyn. Astr. 90 , n. 3-4, 245-266 (2004) 

[4] A. Celletti, G. Pucacco, D. Stella, Lissajous and Halo orbits in the restricted three-body problem, 
Preprint (2014) 

[5] C.N. Douskos, V.V. Markellos, Out-of-plane equilibrium points in the restricted three-body problem 
with oblateness (Research Note), Astronomy and Astrophysics 446 , 357-360 (2006) 

[6] C. Froeschle, E. Lega, R. Gonczi, Fast Lyapunov indicators. Application to asteroidal motion, Cel. 
Mech. Dyn. Astr. 67 , 41-62 (1997) 

[7] H. Goldstein, Classical Mechanics, Addison-Wesley, 3rd edition (2001) 

























QUALITATIVE AND ANALYTICAL RESULTS OF BIFURCATIONS TO HALO ORBITS 


29 


[8] G. Gomez, J.M. Mondelo, The dynamics around the collinear equilibrium points of the RTBP, 
Physica D 157, 283-321 (2001) 

[9] A. Jorba, A. Farres, Dynamics of a solar sail near Halo orbits, Acta Astronautica 67, n. 7-8, 979-990 

( 2010 ) 

[10] A. Jorba, A. Farres, On the high order approximation of the centre manifold for ODEs, DCDS-B 
14, n. 3, 977 1000 (2010) 

[11] A. Jorba, J. Masdemont, Dynamics in the center manifold of the collinear points of the restricted 
three body problem, Physica D, 132, 189-213 (1999) 

[12] J. Laskar, Frequency analysis for multi-dimensional systems. Global dynamics and diffusion, Physica 
D 67, 257-281 (1993) 

[13] J. Laskar, C. Froeschle, A. Celletti, The measure of chaos by the numerical analysis of the funda¬ 
mental frequencies. Application to the standard mapping, Physica D 56, 253-269 (1992) 

[14] C.D. Murray, S.F. Dernrott, Solar system dynamics, Cambridge University Press (1999) 

[15] P. Oberti, A. Vienne, An upgraded theory for Helene, Telesto, and Calypso, Astronomy and Astro¬ 
physics 397, 353-359 (2003) 

[16] A. Marchesiello, G. Pucacco, Relevance of the 1:1 resonance in galactic dynamics, Eur. Phys. J. 
Plus 126, 104 (2011) 

[17] A. Marchesiello, G. Pucacco, Equivariant singularity analysis of the 2:2 resonance, Nonlinearity 27, 
43-66 (2014) 

[18] G. Pucacco, A. Marchesiello, An energy-momentum map for the time-reversal symmetric 1:1 reso¬ 
nance with Z 2 x Z 2 symmetry, Physica D 271, 10-18 (2014) 

[19] R.K. Sharma, P.V. Subba Rao, Collinear equilibria and their characteristic exponents in the re¬ 
stricted three-body problem when the primaries are oblate spheroids, Celestial Mechanics 12, 189- 
201 (1975) 

[20] V. Szebehely, Theory of Orbits, Academic Press, New York and London (1967) 

Departments of Mathematics, University of Roma Tor Vergata, Via della Ricerca 

Scientifica 1, 00133 Roma (Italy) 

E-mail address : kayleigh.s@hotmail.it 

Departments of Mathematics, University of Roma Tor Vergata, Via della Ricerca 

Scientifica 1, 00133 Roma (Italy) 

E-mail address: ceccaron@mat.uniroma2.it 

Departments of Mathematics, University of Roma Tor Vergata, Via della Ricerca 

Scientifica 1, 00133 Roma (Italy) 

E-mail address: celletti@mat.uniroma2.it 

Departments of Physics, University of Roma Tor Vergata, Via della Ricerca Scien¬ 
tifica 1, 00133 Roma (Italy) 

E-mail address: pucacco@roma2.infn.it 



